Driven response of time delay coupled limit cycle oscillators 
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We study the periodic forced response of a system of two limit cycle oscillators that interact 
with each other via a time delayed coupling. Detailed bifurcation diagrams in the parameter 
space of the forcing amplitude and forcing frequency are obtained for various interesting limits 
using numerical and analytical means. In particular, the effects of the coupling strength, the 
natural frequency spread of the two oscillators and the time delay parameter on the size and 
nature of the entrainment domain are delineated. For an appropriate choice of time delay, 
synchronization can occur with infinitesimal forcing amplitudes even at ofi'-resonant driving. The 
system is also found to display a nonlinear response on certain critical contours in the space of the 
coupling strength and time delay. Numerical simulations with a large number of coupled driven 
oscillators display similar behavior. Time delay offers a novel tuning knob for controlling the 
system response over a wide range of frequencies and this may have important practical applications. 



PACS numbers: 05.45.+b, 87.10.+e, 44.27 
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I. INTRODUCTION 

Coupled limit cycle oscillators have been extensively studied in recent times to understand synchronization 
phenomena in various physical^ chemical and biological systems U i 13. | 3t \ ^ l a , |g . Examples include coupled lasers 
3,1 1 , , coupled magiietroiis Fril. arrays of Josephson junctions |lll Il2l Il4l llfil 11^ . coupled chemical reactors 
18l Il9l I20I I21I I23. l25j|. and coupled arrays of biological cells |26l l27l l28|. Many of these practical systems 

are also often subject to external driving forces. For example, the day and night variation of the solar radiation 
input provides a natural periodic forcing of many biological and ecological systems. Electronic pace-maker devices 
implanted in the human body for regulating cardiac rhythm, driven chemical oscillations in industrial reactors, 
and control of chaos through periodic forcing of mechanical systems are a few other common examples. Another 
ubiquitous feature of natural systems is the presence of time delay in the mutual interaction between their component 
elements 293 ,30J- Such delays are normally associated with finite propagation times of signals, finite reaction times 
in chemical systems, or individual neuron firing times in neural networks. Time delay introduces interesting new 
features in the collective dynamics of coupled limi t cy cle oscillators as has been pointed out in a number of past 
studies [3l|33,|3l0IM|3i|3l|3l|3i|4fl|4l|4l|H 

The characterization of the driven response of coupled oscillator systems has important practical applications, and 
has been carried out in the past for a number of simple systems pTI . lisj . A well known case is the driven Kuramoto 
model where the collective phase evolution of a group of weakly coupled limit cycle oscillators due to periodic forcing 
has been examined [4^ . However for many practical applications the approximation of weak coupling is not valid. 
For such cases where the coupling is strong, one needs to take into account the amplitude evolution in addition to the 
phase dynamics. Examples include the Hodgkin-Huxley family of equations used in many biological applications and 
coupled Stuart-Landau oscillators which have been investigated in the past in the context of spontaneous collective 
synchronization phenomena JSO,!^. Recently we have generalized the Stuart-Landau model to include time delay 
effects in the coupling terms |35l l3a | . In this paper we study the driven periodic response of such a generalized system 
when it is subjected to an external oscillating force. For simplicity we restrict ourselves mostly to the investigation 
of just two delay coupled oscillators which are each subjected to an identical periodic force. However as we briefly 
demonstrate toward the end of the paper the results hold for a system of large number of oscillators as well. Our 
main objective is to examine the nature of the synchronized response of the system in various parametric regimes 
in the space of the coupling strength, natural frequency spread, time delay, external forcing strength and external 
forcing frequency. Using both numerical and analytical methods (wherever possible) we obtain detailed bifurcation 
diagrams to mark the regions of stability of driven synchronized solutions. We also examine the nature of the 
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forced amplitude response and its frequency dependence. The system is found to display a nonlinear response on 
certain critical contours in the space of the coupling strength and time delay. Time delay plays an important role 
in determining the driven frequency response and appears to offer a novel tuning knob for controlling the system 
response over a wide range of frequencies. 

The paper is organized as follows. In Section ^ we write down the model equations for the time delay coupled 
driven system and also recapitulate briefly the salient collective properties of the non-driven, no-delay system that 
was studied in detail by Aronson, et al. [50.] . In this section we also define the particular kind of driven synchronized 
response of the coupled system that is the object of our study in the presence of delayed coupling. We further 
derive the eigenvalue equation to find the linear stability of these solutions. In Section UTTI we first study these driven 
synchronized solutions in the absence of time delay. Analytical curves describing the boundaries of the synchronized 
solutions are presented. We find that when the synchronized state loses its stability, the average frequency < 9 > 
could either grow continuously or make a finite jump depending on the driving strength. The scaling behavior of the 
synchronized responses as a function of F is studied for a special case of resonant driving {Co = Q). Finite dispersion 
among the oscillators' frequencies facilitates the synchronization as shown by the full bifurcation diagrams. In Section 
II VI the driven synchronization mechanism is studied in the presence of time delay. Analytic formulae of the boundaries 
of the synchronized state in (t, K) space are obtained. The nonlinear nature of the driven response on these curves 
is elucidated and its possible implications discussed. For finite dispersion the corresponding bifurcation diagrams, 
and frequency jumps across the boundaries of the stable and unstable synchronization regions are presented. As a 
function of r, the frequency jumps are found to scale quadratically. In Section we briefly present numerical results 
on driven synchronization of a large number of coupled oscillators. Section IVII summarizes our results and makes 
some concluding remarks about applications and future directions of research. 



II. MODEL EQUATIONS 

We begin by writing down the model equations for our driven coupled system, 

Zi = (1 + tu, - \Z,f)Z, + K[Z2{t - r) - Zi] + Fe'^\ (1) 

Z2 = (1 + - \Z2\^)Z2 + K[Zi{t - r) - Z2] + Fe'"*, (2) 

where Zj{— Xj + iyj — rjc'^^) is the complex evolution variable, cui, and lu2 are the intrinsic frequencies of the 
oscillators, K{> 0) is the mutual coupling strength, t(> 0) is the amount of time delay in the coupling mechanism, 
and both oscillators are acted upon by a uniform periodic force, Fe*^*. Let us define two useful frequency 
parameters: u) — {loi +lu2)/2, the mean frequency, and A = — W2I the natural frequency mismatch (spread) 
between the two oscillators. The above coupled system without any external driving (i.e. for F = 0) has been 
the subject of detailed studies in the past. We briefly recapitulate the salient collective features of the undriven 
system. For r = 0, as demonstrated by Aronson et al. |50| the system essentially shows three distinct kinds of 
behavior: (i) phase-locked (synchronized) limit cycle solutions with both oscillators oscillating at lu ^ Cu, amplitudes 
rf = r2 = 1 — K + \/K'^ — A2/4, the phases 9i = ivt — a/2, 62 = tot + a/2 with a — sin^^(A/2iir), (ii) asynchronous 
(phase drift) behavior where each oscillator behaves nearly independent of the other and maintains its own natural 
frequency, (iii) collapse of the limit cycle amplitude to zero (amplitude death). A composite phase diagram in the 
K — A space delineating the parametric regimes for the occurrence of these collective states has been provided by 
Aronson et al. |5Clj| and we reproduce their diagram as Fig. ^ here. The phase drift behavior (region I) occurs for 
all values of A > as long as the coupling strength is weak: K < mm{2K, 1}. For small values of A phase locked 
behavior appears after the coupling strength crosses a threshold value (region II) . When both the frequency disparity 
(A) and the coupling strengths are large (i.e. in the region 1 < K < (1 + A^/4)/2, and A > 2) amplitude death occurs 
(region III). The transition from amplitude death to phase locked behavior across the boundary of regions III and II 
is in the nature of a Hopf bifurcation. In the presence of time delay, as earlier shown by us [351 136| , some significant 
modifications occur in the collective dynamics of the system and in the nature of the phase diagram. For example, 
the collective frequency of the phase locked state is no longer at w but at a frequency which has a dependence on 
the time delay parameter r. In fact this dependence is in the form of a transcendental relation giving rise to the 
existence of multiple frequency locked states. Another significant change is in the topology of the amplitude death 
region which can now extend down to the A = axis for some values of r and also show multiple-connectedness. In 
the presence of time delay, identical oscillators can phase lock to either in-phase or anti-phase locked states. Some of 



3 



8 











/ 

/ 

/ 


T 
1 


III 

A 
r\ 


/ III 

/ ^."-^ / 




/ 

/ 

/ 

/ ^ 

/ 

/ yX 

/ 
/y 


iS: = (l+A2/4)/2 




1 


II 

1 1 1 







K 



FIG. 1; Bifurcation diagram of Eqs. 11I2II for r = 0, and F = 0. 



these delay induced properties have also been verified experimentally |52l |5; 



We now wish to examine the driven response of the coupled system discussed above. Note that for finite F the 
trivial state (Z\ = Z2 — 0) is no more a solution of Eqs. (|1I2|I . Thus the state of amplitude death no longer exists. 
However the parametric region corresponding to the death state can now sustain oscillations that are externally driven 
and the nature of this driven response and the domain of stability of these oscillations are the object of our study. 

Our study is focused on a particular class of driven solutions in which each of the coupled oscillators oscillates with 
the external frequency. In other words the oscillators are not only synchronized with each other but also with the 
external frequency. The synchronization is restricted only to the frequency and there can be a finite phase difference 
between the oscillations of the two oscillators as well as with the phase of the driving force. These particular solutions 
can therefore be expressed in the form, 



(3) 



and 9i^2 = + ai^2, where i?i,2, and ai^2 are all real constants. 

Substituting these solutions Q in and JSJ, and separating the real and imaginary parts, we arrive at four 
transcendental equations that determine the amplitudes -R1.2 and the phases ai,2 for a given set of F, il, and r: 



{1- K - Rl)Ri^ 
{oji - - 

{l-K~Rl)R2^ 
{lU2 - n)R2 



KR2 cos + F cos ai — 0, 

- KR2 sin — F sin ai = 0, 
KRi cos A_ + F cos a2 = 0, 

- KRi sin A- — F sin 02 — 0. 



(4) 
(5) 
(6) 

(7) 



where A± = —fir ± (a2 — ai). Notice that the frequency of the system is specified by the external driving frequency. 
The above set of equations can possess a large number of roots (oscillatory states) particularly when r is finite. Except 
for simple limiting cases (e.g. r = 0), these solutions can only be obtained numerically. In addition we also need to 
know the region of stability of these solutions in the (A, K) plane. Following the standard procedure of linearizing 
around the equilibrium solutions, it is easy to write down the stability matrix for the above solutions as. 



M 



Pl + ILOl 

BC 




where pi_2 = 1 — K—2R1 2, wi,2 
equation is then given by. 



Ji.2-n, B = iCe"-^^, C 



—Ri'Ei 

Pl - i<^i 


B/C 

At 



BC 


P2 + 

— i?2 / E2 





B/C 



, and El 2 



— i?2 E2 
P2 - . 

i2ai,2 



(8) 



The corresponding eigenvalue 



L"* + asL^ + a2L^ + aiL + oq = 0, 



(9) 
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where A is the eigenvalue, L = l-fC-A, 03 = -4i?l^-4i^2^ 02 = -^^GjI^^?, Ri"^ + GjI-C'^B'^ + ?, R2'^ + l&Ri^ , 
ai = 2 C^B^R^" - 12 Ei^iJs' + 2 C^^^^i^ - l^^^ + 2 - 4 UIR2' - - 4 R.^Cdj + iC^B^[o2 + 2 + 

iC'^B'^Oji - 12i^l'^i?2^ and ao = cZ'Jwl - AC^B^Ri^R2^ + B^ - 2iC^B^CoiR2^ - 2iC'^B'^ Ri^Cj2 + 2'-^^^^ + 

gnp^ + 2 i^^^^J^ ~ 4 - E2R,-^B-R^ ^ 3 ^^4-2 ^ 3 ~ 2^^4 ^ ^2^2^^ ~^ _ ^^^'^^^^^^^ + 9 R^^ R2\ 

In the following sections we will solve the set of equations Q-CJ in different limits to obtain various synchronized 
states and use equation to determine their domains of stability along with a discussion on other possible solutions 
admitted by the system. 



III. NATURE OF SYNCHRONIZATION IN THE ABSENCE OF TIME DELAY 



In this section we investigate the effect of force on the coupled system in the absence of time delay. We start our 
analysis by considering the case of identical oscillators. When the oscillators are identical, i.e. (wi — uj2 = w), and 
undriven, i.e. F ~ Q, the system admits two kinds of solutions: (i) in-phase solutions characterized by Z\ = Z2 = Z ^ 
and (ii) anti-phase solutions characterized by Zi = — Z2 = Z. We will study here the in-phase solutions under finite 
force. Under the assumption of the in-phase solutions, the study of Eqs. JQ) and (|2Jl becomes identical to studying 
a single oscillator with a delayed linear feedback and constant external forcing. Thus, it is sufficient to focus on the 
following simple equation: 



Z(t) = (1 + - \Z(i)\^)Z{t) + Fa 



'.nt 



(10) 



or under the transformation Y{t) — Z{t)e~ 



Y{t) = {l + iu-\Y{tf)Y{t) + F, 



ill) 



where uj = uj — Such a set of equations in the context a single relaxation oscillator under an external force was 
earlier studied by Guckenheimer and Holmes p^ . and as shown by them, the system admits two kinds of solutions: 
one a stationary state of the above equation that corresponds to a periodic solution of Eq. (|10|) . and the other, a 
periodic solution that corresponds to a solution of Eq. 11U|) with two frequencies. Both these solutions in our present 
context would correspond to the in-phase oscillations of the two identical oscillators represented by Eqs. pi2l) . The 
stable fixed point of Eq. would correspond to a synchronized state of the oscillators that is also synchronized 
with the external frequency, whereas the periodic solution of Eq. Hll|l would represent a synchronized state of the 
oscillators that is not synchronized with the external frequency. 

The stationary solution of Eq. (|ll|l is the synchronized solution of the system. Substituting Y = i?e'" in Eq. 111() 
and separating the real and imaginary parts one arrives at the following two equations: {1 ~ R^)R + F cos a = 0, and 
CjR — Fsina — 0, which can be algebraically simplified to arrive at the following expressions for determining R and 



R"" ~ 2i?* + (1 + Gj')R^ 
sin-i ^ if 
TT - sin-i ^ if 



'F^ = 0, 
l<i?2, 
1 > R\ 



(12) 
(13) 



Eq. H12() produces a curve R = j{u}, K, fir, F) which can be multiple-valued. The number of real roots for R^ now 
ranges from one to three. This number is decided by the following factor: 



G 



1 

27 



B-^ - B^ 



9F^B 



27 



where B = 1 + Q'^. There will be one real solution if G > 0, three real solutions either if G = (in which case at 
least two of them are identical) or if G < 0. In particular, for uj — and r = 0, there will be one real solution when 
F > two solutions at the points F = and F = and three solutions when < F < 

3v3 3v3 oVo 

The stability of these fixed points must be determined starting from matrix |(SJ) by inserting uji =0^2 = lu. This in 
the present case simplifies to the following four equations: 



-.21 1/2 



A = l-2i?2±{i?4-w2} 
\ = 1-2K -2R^ ±{R'^ - w^} 



2-1 1/2 



(14) 
(15) 
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FIG. 2; Bifurcation diagram of the in-phase solutions of identical oscillators with no delay. The shaded region represents the 
stability region of Zi.a = i?e''"'+"\ 

We next examine these equations to determine the regions of stability of the fixed points. The curves that define the 
stability region of the synchronized state are similar in form to those discussed by Guckenheimer and Holmes |47| 
in their study of a single driven Hopf bifurcation oscillator. However, we would now have an additional set of three 
curves arising due to the higher dimensionality of the problem. (All the critical curves Ti, . . ., Tg are derived in the 
Appendix.) In Fig. |21all these critical curves are plotted for two different values of K. Notice that for a) = 0, the 
eigenvalues become {1 — R^,l — 3i?^}. So the branch of the fixed points (7) falling above i?i_2 = 1 is stable for all F, 
and that below is unstable. This lower boundary of stability is given by The two unstable branches arising due 
to multiplicity of 7 merge on and disappear above Fi. At large values of Co, the lower boundary of stability region is 
F3. When 1/V3 < |w| < 1/2, the curve 7 acquires stability in its multi-valued region leading to bistability of fixed 
points. The shaded region in Fig. |2] thus represents the stability region in which the system has at least one stable 
fixed point. 

The stability analysis also produced an additional set of three curves F4, F5 and Fg which, for K = merge with 
the curves Fi, F2 and F3 respectively, and do not exist for K > 0.5. These indicate the intersections or mergers of 
other possible solutions of the system with the in-phase solutions studied above. A closed form for these solutions is 
not possible. But all such existing solutions can, in principle, be determined from the algebraic relations Zi,2 = by 
setting 071 = C1J2 = 

A. Amplitude response and frequency jumps 

We now study the response of the system as a function of F. The response of the system in the region of stable 
synchronized solution is given by Eq. (|12|l . In particular notice that since the basic oscillator has the limit cycle 
solution with a fixed amplitude of unity, the response is always linear for small amplitudes as a function of F. And 
at large amplitudes, the response is nonlinear: R oc F^^^. The in-phase oscillations will be locked to the external 
frequency with a phase difference a = when lj = (i.e. when they are driven resonantly), and with a phase difference 
a = TT on the curve F = lu. 

Since the stable region for small F is accessible only for — fij = 0, it is only when the system is driven with uj 
that a synchronization with the external frequency occurs for small F. Otherwise the system always goes to a two 
frequency state. 

It is now left to determine the nature of the system in the region where the symmetric in-phase solution is unstable. 
Of particular interest is the average frequency of the system when compared to the frequency of the driving. In Fig.|21 
the quantity < 9 > /il is plotted as a function of \u! — il\ while keeping fl = 10. The plot shows that this quantity has 
a finite jump across F3 and no jump across F2 which is consistent with the fact that on the critical curve F3 there is a 
finite imaginary quantity of the eigenvalue, where as on F2, the imaginary quantity is zero (and hence the frequency 
is equal to the external frequency). In fact the jump (J) is directly related to the imaginary quantity on the critical 
curves across which the transition takes place: 
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FIG. 3; Frequency transitions for different value of F as \uj — fl\ is increased. 

Finally it should be remarked that for identical oscillators, region I (see Fig. ^| corresponding to incoherent states 
(for F = 0) does not exist. Hence the above transition curve is the only relevant one. 

B. Effect of finite dispersion 

Finite dispersion facilitates the synchronization mechanism as we describe in this section using bifurcation diagrams. 
As seen in Fig. ^ regions A and B together make up the amplitude death region in the absence of an external force. 
The region is bounded by the curves K — 1 and K = 7(A) = ^(1 + A^/4). The eigenvalues that determine the 
stability of this steady state are given byA = 1 — K ± iu) ± \J K'^ — A^/4. The regions A and B are separated by 
the curve K = A/2. Note that the system in the death state shows only one collective frequency f = ui, which 
occurs in region B. In region A, i.e, when K < A/2, the oscillators spiral to the death state with independent 
frequencies: /i^2 = w ± \/ A'^/A — K"^. In the presence of the external frequency interacts with /i_2 and / in the 
regions A and B respectively. Since the system offers damping to these modes, the system as a whole acts like a 
damped nonlinear oscillator with these characteristic frequencies. Thus the frequency of oscillation under forcing in 
both the regions is that of the external driver. The region of stability of this driven synchronization (i.e. the solution 
Zi^2 = i?i,2e'''^*''""^'^'') will have to be self consistently determined using the eigenvalue equation described in the 
previous section. Hence in the presence of F , the region of death now transforms into an oscillatory region in which 
both the oscillators synchronize with one another and also synchronize with the external frequency. The actual region 
of stability of the synchronized solution envelops the earlier death region. Before we present the stability region of this 
solution, we make some observations about the response of the individual oscillators. Without loss of generality, let 
us assume uji > uj2- If > w, then Ri > i?2- If 11 = u), then Ri = i?2- And if O < w, then Ri < R2. At the resonant 
driving, the sum of the phases of the oscillators with respect to the external driver will be 0. i.e. ai + 0:2 = 0. 

We now briefly describe the solutions when the oscillators are driven resonantly, i.e. when Q — uj. At resonant 
driving, it is possible to arrive at a simpler set of algebraic equations to determine the amplitudes and the phase of 
this oscillatory state. Substituting Ri = R2 = R and a2 — —ai in Eqs. (|4I7|I . the following two relations can be 
written down that describe the amplitude and the phase of the system. 

R = F sinai/[LJi ~ i^T sin(2ai)], 

where ai is determined from the following functional relation: 

g{ai) = F^ sin^ ai - {uji - ii'sin(2Q!i))^[l - K + ii'cos(2Q!i)] sinai - {uji ~ ivT sin(2ai))^ cosai = 0. 

The response grows linearly with force. The dependence of a however is not apparent from the equation. Notice that 
the resonant point does not necessarily fall on the critical point of Hopf-bifurcation. The response of the oscillators 
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FIG. 4: Dependence of R and ai on K and F. (a) and (b) are plotted for F — 1, and (c) and (d) are plotted for K = 5. The 
other parameters for all the plots are = 10, A = 18, SI = 10, and t — 0. 




and the phases are shown in Fig. ^ as a function of both K and F. For small F, R grows linearly, and for large F, 
it grows sub- linearly according to i? cx F^/'^. The scaling of the amplitude as a function of F in a linear fashion is in 
accordance with the fact that inside the region of death state, the driven oscillator shows a linear response. 

The coupled system shows an interesting response as a function of K. Increasing K will sweep the stable syn- 
chronized region. For small K , the system shows incoherence, which is marked with dashed lines. The synchronized 
response loses stability at large values if in a Hopf bifurcation. Close to this critical boundary, the phase of of the 
first oscillator as shown in the figure damps as 1/F. However, the amplitude has a very nonlinear behavior as a 
function of F. It can be considered as linear only in the middle of the stability region. Note, however, that since 
the region encompasses a Hopf bifurcation curve, 7(A), the response of each of the oscillators becomes nonlinear and 
proportional to F^/'^ on this curve. This characteristic nonlinear response at the Hopf bifurcation threshold has been 
invoked earlier to model the behavior of inner hair cells in human cochlea |54l l55j | . In the present context of two 
coupled oscillators, we have numerically verified this feature. The phase ai, which is a measure of the phase difference 
the first oscillator maintains with the driver, is nonzero, and is nearly constant for all weak forces, but damps with 
an approximate scaling of F~^ '^ for large values of F . The dotted curve drawn is only an indication of the power 
law. This behavior was not noted in literature before and may have some practical applications. The coefficients of 
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(a) 




FIG. 6: The time averaged quantities < ^1,2 > and < ri,2 > plotted in (a) and (b) respectively. 

the eigenvalue equation ^ reveal that the stability is a function of the product term Qt. Thus the stability region 
in the absence of time delay is independent of fl and only depends on F. We now present the stability regions of 
this synchronized state in the parametric space of {K,A) for different values of F in Fig. [SJa). For weak forcing, 
the oscillators can be moved out of the synchronized regime by strongly coupling them. But as the forcing strength 
increases, much stronger or much weaker coupling is required to destabilize the synchronized state. At weaker cou- 
plings, it is the dispersion (A) that is responsible for instability, and at stronger couplings, the system is more likely 
to get locked to its average frequency than to the external frequency. The relation between F and \uj — n\ is shown 
in Fig. [Sfb). At larger values of K, a strong force is required to synchronize the system to the external frequency. 
When these solutions become unstable, there are two other possible solutions that are similar in character to the 
synchronized, and the incoherent states of the non-driven system. For any of the parameters corresponding to the 
bifurcation diagrams in Fig.[SJa), we note that the system's average frequency is locked to the average frequency Hi 
of the oscillators for small value of A. A numerical plot of the averaged frequency is shown in Fig.lH^a) which reveals 
that the transition from one coherent state to another is discontinuous. The nature of bifurcation at this transition 
is supercritical Hopf as can be seen from Fig. |HIb). 

IV. NATURE OF SYNCHRONIZATION IN THE PRESENCE OF TIME DELAY 



In this section we investigate effect of time delay on the forced synchronization of the coupled system. We again 
start our analysis by considering the instructive case of identical oscillators. We will study here the in-phase solutions 
under finite force. As in the previous section, under the assumption of the in-phase solutions, the study of Eqs. 
and |5J becomes identical to studying a single oscillator with a delayed linear feedback and constant external forcing. 
Thus, it is sufficient to focus on the following equation: 



Z{t) = {I- K + iuj- \Z{t)f)Z{t) + KZ{t - t) -f Fe 
or under the transformation Y{t) = Z(i)e~*^*, 

Y{t) ^{l-K + iCo- \Y{t)\^)Y{t) + Kc-'^^^Y{t - t) - 



(16) 



(17) 



where uj = uj — fl. Substituting Y = i?e'" in Eq. (|17|l and separating the real and imaginary parts one arrives at the 
following two equations: 

{I- K - R'^)R + KR cos(rJr) + F cos a = 0, 
CjR - KR sin(f^T) - F sin a = 0, 

which can be algebraically simplified to arrive at the following expressions for determining R and a: 

aiR^ - = 0, 
if c < 



F 



R'' + 02 i?* 
TT — sin 



if 



c> R\ 



(18) 
(19) 



where b 



ifsin(riT), fli = b'^ 



0-2 



-2c, and c = 1 — K + Kcos{ilT). Eq. IjlSfl produces a curve 



R = j{uj, K,flT, F) which can be multiple-valued. The number of real roots for R^ now ranges from one to three. 
The stability of these solutions that can be determined from M are governed by the following eigenvalue equations, 



A = 1 - X - 2i?^ 



SkKe 



-Xt 



(20) 
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where Sk = ±1 and s = ±1. These two signs must be taken in all the permutations. Thus we have four transcendental 
equations. It is possible to treat all these equations together. Just as we did in the no-delay case, we consider here 
also two cases (i) > and (ii) i?^ < In the first case, in order to arrive at the critical curves across which 
transitions of eigenvalues take place, we set A = if3. This gives the following two equations: 



1~K - 2R*^ + sVr*^ ~Cb^ ^ ~SkK coaiPr), (21) 

(3 = SkKsiniPT). (22) 



Note that /3 = is a solution of the Eq. (|22|l for any value of K and t. This branch corresponds to the critical 
curves that existed for the case of t = when both s and Sk are considered independently. For a second branch of 
solution for Eq. (|22|l to exist, the necessary condition is Kt > 1. Again from Eq. (|22|l . it can further be concluded 
that the critical values of Kt beyond which the second branch exists are determined from the zeros of the function 
g = cos{y/ K'^t'^ — 1) + 1/ [skKt). A numerical analysis of the equation g = gives the critical values: Kt = 4.6034, 
if — 1, and Kt = 1, if = — 1. A set of critical curves for the second and higher branches are obtained by 
explicitly considering s = ±1 of the newly generated branches. Also note that the response now is a function of r. 
Hence the stability can be affected by the /3 = branch itself. An equation for R^ can be written down by eliminating 
fi as follows: 



P^lK^- 



l-K-R*"^ + s 



■[r*'-q'} 



1/2- 



1/2 

> , (23) 



R*^ 



^(^1- K + Sfeif cos(/3r) + sVR*^ - (D^) . (24) 



The critical curves in {F, ui) plane cannot in general be written down in a closed form as we did for the case of r = 0. 
The conditions that exist on K, and uj must now be determined using the two transcendental equations above. These 
equations can be used along with Eq. ((T51) to numerically plot the critical curves in (F, uj) plane by setting i?^ = R*^. 
We can easily identify that the curves that are counterparts of Fi and F2 of the no-delay case are generated here by 
setting Sfe = 1. The critical curves corresponding to /? = can, however, be determined in closed form. Following a 
similar analysis as outlined in the Appendix, these curves are obtained as 

r'l = If = f{p+)\p+ = ^ + lVl- 3CuA,u;' e [0, i], (25) 



3 3 



= = f{p-)\p+ = 1- JVT^s;^} e [0, i], 



(26) 



where f{p) = {p^ + a2p^ + aip}^^^. These are the only critical curves that exist under this case as long as Kt < 1. 
If this condition is violated pairs of sets of such critical curves exist. 

In the second case, i.e. when R^ < \uj\, on the critical curves, we again set X = i/3 to arrive at the following set of 
equations: 

1-K- 2i?*^ = ~SkK cos(/3r), (27) 
f3 + SkK sin(/3T) = s^Jr*'^ -Oo^. (28) 

Note that since the eigenvalues occur in complex conjugate pairs, the second of the above two equations yields the 
same set of curves for both s = ±1. So there will be two classes of curves corresponding to s^: = ±1. The above two 
equations can further be simplified to yield 



/? = /3± = ^Gj^ - R*^ ±^K^ - (1- K - 2i?*2), (29) 
R*^ ^^(1- K + SkK cos(/3t)) . (30) 

The second set of critical curves, F3, in {uj, F) plane are obtained using the above two equations along with Eq. H12I) by 
setting R^ — R*^ . We can again identify that the counterpart of the curve F3 of the no-delay case is obtained here by 
setting Sfc = 1. As r is increased, the parameter region falling below the curve F3 becomes more and more susceptible 
for synchronization, and in fact for a range of value of r the curve makes intersections with _F = line as shown in 
Fig. [3 The region marked I is the region of stability of the synchronized solution. At the region of intersection of 
the three curves, there is also a possibility of more than one synchronized solution coexisting. At higher values of 
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x=0.15, D.=l, K=20 




FIG. 7: The critical curves F'l, r'2, and T'^, in the presence of r at A" = 20, fl = 1, and r = 0.15. Region I represents the 
synchronized region. 

time delay, a stronger force is required to make the oscillators synchronize with the external driver. The other set of 
curves that are generated using the second sign = — 1 represent intersections of the response of the system with 
other unstable solutions the system possesses. We do not discuss them here. This figure depicts one of the important 
results of our present paper. In the presence of appropriate time delay the region of forced synchronization can be 
extended to continuously connect with F = axis. This raises the interesting possibility of achieving synchronization 
far from resonant driving with a minimal (nearly infinitesimal) force. 

The role of time delay can be more clearly appreciated by determining the required critical curves in the plane of 
(t, K) . By following the usual analysis of such characteristic equations and properly choosing the correct signs for the 
cosine function above, the critical curves in the plane of {K, r) can easily be derived. When R*^ > \uj\, the equations 
for R and f3 at criticality, Eqs. (|^ and can be used to derive an equation for r in terms of K: 



rnr — cos 



K 



(31) 



where A± ^ 1 ~ K - 2R*^ ± \/r*^ — uj^. On both of the above two curves, t±, the nature of the transition of the 
eigenvalues is given by: 



^^ReiX) 



Re{\)=0 



> 0. 



(32) 



Thus these two curves in the (r, K) space represent critical curves across which a pair of eigenvalues makes transition 
to acquire positive real parts. 

In the second case, i.e. when R*^ < \uj\, the valid equations are Eqs. (|27|l and H28|l . By again following the 
standard methods to derive the critical curves, we arrive at the following set of two marginal stability curves: 



rnr — cos 



T2 



-1 { l-K-2R'^ \ 



R* 



nil + cos 



K'^- {l-K -2R*'^Y 

-1 ( 1-K-2R' 



- R*^ - Jk^ -{l-K- 2i?*2)2 



(33) 



(34) 
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T 

FIG. 8: Jumps of the average frequency across the critical boundary separating synchronous and non-synchronous regions. 
K — 0.4, F = 0.6, Q. — 10. Bullets are numerically found values. The curve is a best fit with S.eSSr^ + 0.043. 



The nature of transition of the eigenvalues across these critical curves is given by: 

-^i?e(A) pOonn, 
dr R.e{X)=o l<0 on T2. 

The curves t± and ti_2 completely describe the stability transition curves in (r, K) space for a given set of F and 
to. The ordering of these curves depending on F and uj will eventually determine the curves that enclose the stable 
regions. 



A. Amplitude response and frequency jumps 

We now examine the response of the system. Notice from Eq. H18|l that the oscillators respond linearly with F 
for small amplitudes (i.e. when the F^ term is balanced by the aiR^ term). For larger amplitudes the response 
becomes nonlinear as the balance is provided by other terms. However the response can be highly nonlinear even for 
small amplitudes on certain contours of the parametric space where 02 = 0, and oi = 0. Under these conditions the 
following two equations are obtained (provided that Qt 7^ 2mr): 

cos(17t) = (1 - K)/K, 

and 

sin(f7r) = {uj~fl)/K. 

Using the above two equations we can derive the conditions relating K and Qt for a given uj such that the response 
of the system for small amplitudes has a nonlinear character {R oc F^/"^): 

K = (l+cZ'2)/2, (36) 

Qt = 27i7r + cos"i (^Y~r^) ' ^^"^^ 

In the [F, Co) plane, the corresponding contours can be derived from Eq. (|18|l by making use of the above two relations. 
We will discuss this phenomenon in more detail later. At the boundary of stability region, the oscillators make a 
transition in their average frequency from their intrinsic average frequency to that of the driver. This discontinuity is 
a function of r and is numerically found to be quadratic in its dependence on r as shown in Fig. |H1 



12 



a =10, (0=10.A=18, K=15, F=1 




o.Sir/si It/a i.5it/a 2it/a o o.5it/a i/a i.5ir/a 2it/a 

T T 



FIG. 9: Amplitudes and phases of the osciUators as a function of r. 

B. Effect of finite dispersion 

We now study the effect of finite time delay and the external force on non-identical oscillators. Our primary objective 
is to examine the stability properties of the synchronized state Z\^2 — i?i,2e**-^~''"^-^-' which is also synchronized with 
the external frequency. We use Eqs. H4I7(I to determine the responses i?i,2, and ai_2 and then find the eigenvalue 
spectrum to determine their stability. The response due to time delay shows considerable deviation from the no-delay 
case. Under certain conditions, both the oscillators could have identical amplitudes. Similar to our studies in the 
previous section, we assume that Ri = R2 = R and ai -I- a2 = 0. This ansatz leads to the condition that Or = tt. 
Substituting these relations in Eqs. (j4l7|l . we arrive at the following two simplified equations for the responses: 

R 

cosai = ~—{l-K-R^-Kcos{2ai)), 
R 

sinai = — (cji + i^T sin(2a)) , 
F 

which can further be simplified to arrive at an equation for i? as a function of ai as 

R^ Fsinai/[Cji+ Ksin{2ai)\, (38) 

and a functional relation to determine ai as 

h{ai) = F'^ sin^ai - (wi -I- iv:sin(2ai))^[l - K - if cos(2ai)] sinai - {Cui + ii'sin(2ai))^ cosai = 0. 

The above expression for cos ai again indicates that for this special case the response of the oscillators is nonlinear 
according to i? oc F^/'^ at large values of F. 

However, the amplitudes of the oscillators are in general not identical. The response of the oscillators in the 
synchronized state is unique only when < ilr < 27r. In Fig. El the amplitudes and phases are shown as a function 
of r. For a given r, stronger coupling among oscillators decreases the amplitudes. At Jlr = tt, the oscillators are 
locked with the external driver with a maximum phase difference. For a given set of parameters the external forcing 
will widen the region of stability of the synchronized solution. In Fig. llOf a') the response of one of the oscillators is 
plotted for different values of F. The boundary of the stability regions of the death state in the absence of external 
driving falls inside the stability region of the synchronized state. In Figs. [TUT b). (c), and (d), the stability region of 
the synchronized state is shown in (i^T, A) space for three different values of the driving frequency Q. and a fixed set 
of F = 1.0, T = 0.05, and u) = 10. In each of these figures, the dashed curve indicates the marginal stability curve 
for F = 0. As is evident the stability region of the tuned state is sensitive to Vl and appears to increase in area 
for slow driving. The difference between the effects of force (as seen from Fig. O and that of time delay is evident 
here in the bifurcation diagram. Time delay breaks up the degenerate bifurcation point where as just the presence 
of force alone would not remove the degeneracy. Finally we examine the average frequencies as a function of the 
frequency dispersion. In the absence of time delay, as we noted in the previous section, there is a transition of the 
average frequencies to the external frequency as A is increased. However, in the presence of time delay, the intrinsic 
frequencies are affected by r differently, and thus there could be connected regions of r where the oscillators can be 
found in two different locked states as shown in Fig. ^] Here the first oscillator with small frequency (left curve in 
the frequency plot) is more susceptible to the external driving and gets locked to it at a smaller value of A while the 
second one is still locked to the branch of the average intrinsic frequency. The corresponding amplitudes are shown 
in Fig. IllT b). This feature is absent for r = and may be found in multiple connected regions of r. 
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FIG. 10: (a) Response of the first oscillator for various as a function of K. Stability regions (shaded) of the tuned state in 
{K, A) space for (b) = 12, (c) f2 = 10, (d) f2 = 8, while keeping the other parameters at F = 1, r = 0.05, and a) = 10. 




FIG. 11: Average frequencies and the average amplitudes as a function of A. The first of the two oscillators has smaller 
frequency. In (a) the left vertical curve corresponds to the first oscillator. F = 1, r = 0.2, a) = 10, O = 2, and K = 1.5. 



C. Tuning property of time delay 

Having now found the region of stability of the tuned solutions, we must then examine the response of the oscillators 

in this region and find where the maximum or resonant response occurs and with what functional relation this response 
scales with the driving force, F. The Hopf bifurcation curves of the non-driven system in the presence of finite r have 
a special significance in that on them the response is nonlinear even for small F. This fact was noted earlier in the 
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FIG. 12: (a) Two of the intrinsic frequencies existing in the undriven system. The sohd portions of the curves represent damped 
regions, i.e. the real part of the corresponding eigenvalue has negative sign. The bullets indicate the values of K at which the 
damping of that particular mode is zero, (b) The tuned frequency as a function of r. (c) Resonances on the tuned frequency 
curve, (d) The amplitude response of one of the oscillators on and around the tuned frequency curve. The curve 71 is a straight 
line with slope unity, and 72 is that with 3. 



literature and was used to model the nonlinear compression of sound frequencies by the inner hair cells of the cochlea 
f54, 55). In the presence of r the frequency on such critical curves has a dependence on r. This makes it possible to 
choose a proper time delay such that a nonlinear response at a given external frequency is achieved. 

We explore this interesting feature further by examining the eigenvalue spectrum of the death state in the absence 
of external force. The eigenvalue equation that determines the stability of the origin is 



A=l-A'±y'X2c-2Ar_^2/4±^ (39) 

which can be derived from Eq.l^by inserting R12 = 0. We are now interested in the imaginary parts of the eigenvalues. 
The real parts only reveal the damping of that particular frequency mode when the driving is turned on. A detailed 
study of the marginal stability was earlier carried out in literature. Here we examine the imaginary parts of the 
eigenvalue spectrum (i.e. the intrinsic frequencies of the system) as a function of K. As can be verified from the 
full spectrum, the lower frequencies are least damped and the higher frequencies are heavily damped. So we focus 
our attention on the first two frequencies of the system. Since the eigenvalues occur in complex conjugate pairs, 
it is sufficient to plot the positive frequencies. In Fig. I12f a^. the intrinsic frequencies of the eigenvalue spectrum 
are plotted while showing the nature of the damping existing for each branch. The system responds linearly to the 
external driving in the region of damping. The two frequencies plotted have critical points on them marked a, b and c. 
A nonlinear growth of the amplitude of the oscillators occurs at all the points a, b, and c as a function of the driving 
force at any given value of The higher frequency branch has only one critical point a at which the damping is zero. 
To the left of the point, the frequency is a growing mode and to the right, the frequency is a damped mode. Also note 
that the region left to point a indicates a drift region. The lower branch has two critical points on it: b and c both 
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representing zero growth of the mode. The region left to the point b and that to the right of c are growing modes. 
The critical points a and b are not sensitive to variations in r and thus their frequency span is too limited. However 
the point c is sensitive to r and in fact spans a range of values as r is increased. Thus this point is of interest to us. 
At higher values of A, a large frequency range can be spanned with a short variation in t. In Fig. I12r b) the contour 
of c is plotted in {^,t) space. At this frequency the system responds nonlinearly and can be termed as the tuned 
frequency. In Fig. Ifif c) the frequency response of the system is shown as one traverses along the tuned frequency 
curve shown in Fig. ^Jb) for A = 10. 

Now we show the response of the oscillators at the tuned frequencies. The system responds nonlinearly with the 
amplitude proportional to F^^^. In Fig. I12r d') the amplitude of one of the oscillators is plotted in log-scale for a 
choice of the parameters lying on the tuned frequency curve shown in Fig. ll^f b'): K{— K*) — 16.861411237, r = 0.06, 

— {p,*) — 4.968163, [D = 10, A = 6. The responses around the curve are Hnear, but on the critical curve the 
response is nonlinear and is proportional to F^^'^ . 

V. SYNCHRONIZATION OF A LARGE NUMBER OF COUPLED OSCILLATORS 

We now briefly consider the question of the possibility of synchronizing a large group of globally coupled oscillators 
that have a distribution of frequencies, ^(cli), and are driven by an external periodic force: 

K ^ 

Z,{t) = (1 + ^u^, \Z,{t)\')Z,{t) + -Y^[Z,{t-T)- Z,{t)] + Fe'"*, (40) 

fc=i 

The contribution of the self coupling term that would emanate from the summation above would be ignorable as the 
size of the system becomes large. A full bifurcation picture of these equations is relegated to a future study. Here 
we present the main features exhibited by this system which are common to the case of two oscillators discussed 
so far. g{uj) is assumed to be in the form of uniform density distribution, i.e. g{u;) = l/2j for uj e [—7,7], and 
otherwise. The frequencies are chosen with equal spacing. Since time delay removes the symmetry enjoyed by the 
system, a constant positive value uj is added to the frequency distribution to make all the frequencies positive. A 
convenient macroscopic parameter to characterize the collective properties of a large population of oscillators is the 
order parameter defined by 

1 ^ 

S = R{t)e^m = _J2z,{t). (41) 
i=i 

In Fig. I13r a') the time evolution of = 800 oscillators are plotted in (A, Y) plane in the frame moving with the 
external frequency for F = I, t = 0.18, uj — 1.2, and Q = 1.0. A stationary asymptotic state in the fourth panel 
indicates a complete synchronization of all the oscillators in the population to the driving frequency. Since the forcing 
is strong, the synchronization time is short, and the coupling strength is overcome by F. In Fig. lT^ b) magnitude (R) 
and the frequency ($ of the order parameter are plotted in the moving frame. This confirms the synchrony. 

One of our important results in the earlier sections is the possibility of synchronization for infinitesimal forces in 
the presence of time delay. For this we carry out simulations on A^ = 100 for finite time delay, and finite frequency 
spread for various values of driving force. The results are presented in Fig. ^] as a plot of < $ > /n, vs. r. For a 
strong force {F = I) the horizontal line indicates synchronization. The synchronization occurred here even for t — 
and extends to a very large value of r. As we decrease the force, synchronization at t = is lost, but is attained for 
a range of r values. Even for F = 0.01, this particular choice of 7 shows a range of r values where synchronization is 
found. This confirms our earlier results on two coupled oscillators. Further work on the bifurcation structure of the 
large number of oscillators will be reported elsewhere. 

Finally we look at the frequency transitions from synchronization between the oscillators to that to the driving 
frequency. In Fig. [TsT a) and (b) the average values of the amplitude of the order parameter < R{t) > , and the ratio 
of the average frequency to the driving frequency, < $ > /fi are plotted as a function of the frequency width 7 for 
T = 0, and T = 0.1. At lower values of 7, as is also evident even for two coupled oscillators, the oscillators are not 
synchronized to the driver, but make a transition to the synchronized state at a critical value. At large values of 7 
synchronization loses its stability. The frequency of the collective state at lower 7 is suppressed by time delay. 
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FIG. 13: {&){X,Y) snapshots of time evolution of 800 coupled oscillators under external force, (b) The evolution of the 
amplitude and the derivative of the phase of the order parameter in the moving frame. 7 = 0.6, K = 0.8, F = 1, t = 0.18, 
Q = 1.0, and uj = 1.2. 




FIG. 14: Synchronization regions for different F as a function of r. A'^ = 100, 7 — 0.6, K = 0.8, = 1.0, and lo = 1.2. 



VI. CONCLUSIONS 



We have carried out a detailed theoretical analysis of the driven periodic response of a system of two time delay 
coupled limit cycle oscillators that are subjected to identical external forces. In particular, we have examined the 
stability properties of periodic phase locked oscillations of the system that are synchronized with the external driving 
frequency. Bifurcation diagrams in the parameter space of the coupling strength and frequency spread of the two 
oscillators have been presented and their variation with time delay has been discussed. We have also obtained the 
characteristics of the response function in various limits and highlighted their sensitivity to time delay. We have 
further confirmed that some of the important results of the two oscillator system also emerge in simulations carried 
out on a large number of globally coupled limit cycle oscillators. 

For the sake of completeness and to gain an overall perspective of our results we have also showed the connection 
of our work to earlier investigations of non-driven coupled oscillators, with js^ or without [50] time delay, and to 
driven oscillators without time delay . The parametric region corresponding to the amplitude death of two coupled 
non-driven oscillators has been the primary area of our interest. This is the region that sustains a periodic oscillation 
when the system is driven externally. For our time delayed system the extent of this region and the nature of the 




driven response within it are all sensitive functions of the driving frequency, driving amplitude and the time delay 
parameter. This is evident from the detailed results presented in the various previous sections. In general, the size of 
the stability region of the forced response is larger than the size of the original death region of the non-driven system. 
However the bounding curves of the original region play a special role in that the small amplitude response on them 
is highly nonlinear. 

Two important results of our work are worth highlighting in view of their potential applications. First, by an 
appropriate choice of time delay, it is possible to mitigate the disadvantages of far from resonance driving and to 
achieve synchronization with a minimal of force. Second, for a fixed driving frequency, the driven response of the 
system can be tuned to any desired frequency with a suitable choice of the delay parameter. This property relies on 
the nonlinear response of the system on the Hopf bifurcation boundaries. Since the location of these boundaries (in the 
frequency space) is a sensitive function of time delay it provides for this possibility of using the time delay parameter 
as a control knob for frequency tuning. It would be interesting to look for evidences of these effects operating 
in natural biological and chemical systems or to implement them in appropriately designed electronic receiver systems. 

Finally we would like to remark that in our analysis we have not examined the driven response in the region of 
phase drift solutions. Such a study is presently in progress as is the extension to a larger system of coupled oscillators 
where regions of chaotic solutions could be explored. Future work should also focus on large number of locally coupled 
oscillators. 



APPENDIX: DERIVATION OF CURVES Ti TO Te 



In order to obtain the bifurcation diagram in (F, lu) plane, we must find the critical curves of the above equations 
by setting Re(A) = 0. For the sake of convenience consider two cases: (i) > and (ii) < In either of 
these two cases the real part of the eigenvalues of Eq. I|14(l are always less by a positive value 2K than those of the 
Eq. (|15|l . Thus, at the critical set of parameters where the first two equations have their real parts zero, the real parts 
of the other two equations are negative. Hence, it is obvious that the boundaries of the steady state are completely 
determined by the eigenvalues of Eq. H14|l alone. In fact these two equations determine the stability of the curve 7. 
The other two equations provide critical curves that represent transitions of the eigenvalues. We find that the other 
solutions the system admits will merge with the in-phase solutions at these critical lines. First we derive the critical 
curves under both the cases mentioned above for Eq. (|14() . Then we treat the other two eigenvalue equations. 

Considering Eq. ffT^ with + sign under case (i), Re(A) = 1 — 2i?^ + \/R* ~ ui'^, and Ito(A) = 0. The criticality 
occurs at Re(A) (and thus 1 - 2i?*^ < 0), which results in i?^ = = 2/3 ± ^l/9-a)2/3. Now by making use 
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of Eq. (|12|l . and the condition 1 — 2R*_^ < 0, the following critical curve (Fi) can easily be written down: 



F = 



27 



(l + 9cD') + (1-3.^2)3/2 



1/2 



[i ii 



(A.l) 



Similarly, making use of Eq. (|12|l . and the condition 1 — 27?^^ < 0, the following critical curve (r2) can be written 
down: 



F = 



2 r 

27 



(1 + 9^2) 3-2)3/2 ^ cD2g[0,^] 



1/2 



1, 



(A.2) 



Under the case (i) again, Eq. 11411 with — sign yields Re(A) = 1 — 2R^ — ^/R^~~u^ ^ and Im(A) = 0. The criticality 
occurs at Re(A) = (and thus 1 - 27?*^ > 0), which again results in = R*\ = 2/3 ± ^l/9-a;2/3. Note that 
R = i?!^ fails to obey the condition 1 — 2R > 0. Hence R = i?!^ does not result in any critical curves. But 
R^ — R*_^ results in the following critical curve: 



Ti : F 



2 r 

27 



(l+9c^2)^(^_3~2)3/2 



1/2 



cD^e [0,-]. 



(A.3) 



Under case (ii) Eq. (|14|1 with both signs can be considered together: Re(A) = 1 — 2i?2, and Ir7i(A) = ±\/Cj'^ — i?*. The 
criticality occurs at i?^ = i?*^ = 1/2 with the condition that R^ < \Cj\. This, when substituted in Eq. 112() . results in 
the following critical curve: 



i^^{(^2 + l/4)/2}'/', <i2e[-,oo). 



(A.4) 



Since the imaginary part of the eigenvalue spectrum is nonzero, the nature of the bifurcation across this critical curve 
is determined by 



^ReA 
dF 



ReA=0 W 



-4f r > --d Co' <{, 

i2_l\<0 if Ci^^l^ 



(A.5) 



which indicates that the bifurcation is of the supercritical Hopf kind. 

We now consider Eq. ifTCjl with + sign under case (i): Re(A) = 1 — 2K — 2R' + \/R^ — uj'^, and Im(A) = 0. Let us 
define f{p) = {p^ — 2^^ + (1 + Cj'')p\^^'^, and a — 2K — 1. Following exactly the same method as described above, the 
critical curves can be derived as follows for K < 



r4 : F = {/(p_ 



T,:F^ {f{p+) 



-2a (1)2 _2 , 2 / , 2 /oi 

— - Y y - y}, we [ay A, a73]. 



-2a la'' J;2 ^ „ ^ „ , , 



And using the Eq. I|14|l with — sign, the corresponding curve turns out to be: 

Ti-.F^ {f{p 



—2a la' a)2 



P- = 



g -y}, (i^e [0,aV4]. 



(A.6) 
(A.7) 

(A.8) 



Still under case (i), the nature of the curves under K > ^ differ from the above. The curves are derived using the 
same method as above. Eq. (|14f) with — sign results in 



r,:F = {/(p_) \p_ = ^-^^-^}, Co'e [0, aV3]. 



T,:F = {f{p+) p+ 
and the Eq. H14|l with + sign results in 

r,:F = {f{p+) p+ 



-2a 



2a / a' (1)2 



}, Cj^(.%a'/A]. 



(A.9) 
(A.IO) 

(A.11) 
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Now we are left with case (ii). The equation H14() with both the signs can be treated together: Re(A) = 1 — 2K — 2R^ 
and Im(A) = ±^/R'^ — a)^. On the critical curve, — R*^ = ^ — K. This gives rise to the following critical curve in 
{F,uj) plane: 

r,:F = {i^-K)[^' + {K + ^fy^\ X<i ^2^[i-X,oo). (A.12) 

Since the imaginary part of the eigenvalue is nonzero in this case, the nature of transition of eigenvalues is determined 
by 



dReA 



dF 



-iF f>0 if C^^<\-K-3K^ ^^^^g^ 



ReA=o Cj^ - a- K -3K^)\<0 if Cj^ > \- K -3K^, 



which indicates that the bifurcation involved is once again a supercritical Hopf bifurcation. 
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